################################################################################################
# This script takes the bootstrapped QV data under the rules specified in the QV Welfare  
# Simulations Script and then produces inequality measures and graphs.
# 
# All graphs are exported to SV & QV Recalibration on Pop Imm/QV - Inequality/
# 
# Author: Luis S.
# Last Modified: 7.25.16
################################################################################################

library(dplyr)
library(ggplot2)
library(ineq)
library(stargazer)


setwd("C:/Users/Luis/Dropbox/Research/SV vs QV California Experiment/Experiment Analysis & Results")

QV_IndivUtility_NQV_A <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_NQV_A.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_NQV_B <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_NQV_B.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_NQV_B_CNV <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_NQV_B_CNV.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_NQV_C <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_NQV_C.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_NQV_D_Emp <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_NQV_D_Emp.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_NQV_D_Unif <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_NQV_D_Unif.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_NQV_D_5050 <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_NQV_D_5050.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_NQV_E_Emp <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_NQV_E_Emp.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_NQV_E_Unif <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_NQV_E_Unif.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_QV1_A <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV1_A.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_QV2_A <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV2_A.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_QV1_B <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV1_B.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_QV2_B <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV2_B.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_QV1_B_CNV <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV1_B_CNV.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_QV2_B_CNV <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV2_B_CNV.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_QV1_C <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV1_C.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_QV2_C <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV2_C.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_QV1_D_Emp <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV1_D_Emp.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_QV2_D_Emp <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV2_D_Emp.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_QV1_D_Unif <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV1_D_Unif.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_QV2_D_Unif <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV2_D_Unif.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_QV1_D_5050 <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV1_D_5050.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_QV2_D_5050 <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV2_D_5050.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_QV1_E_Emp <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV1_E_Emp.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_QV2_E_Emp <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV2_E_Emp.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_QV1_E_Unif <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV1_E_Unif.csv", header = T, stringsAsFactors = F)
QV_IndivUtility_QV2_E_Unif <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_IndivUtility_QV2_E_Unif.csv", header = T, stringsAsFactors = F)


QV_RuleAResults <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_RuleAResults.csv", header = T, stringsAsFactors = F)
QV_RuleBResults <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_RuleBResults.csv", header = T, stringsAsFactors = F)
QV_RuleBResults_CNV <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_RuleBResults_CNV.csv", header = T, stringsAsFactors = F)
QV_RuleCResults <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_RuleCResults.csv", header = T, stringsAsFactors = F)
QV_RuleDResults_Unif <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_RuleDResults_UniformVC.csv", header = T, stringsAsFactors = F)
QV_RuleDResults_Emp <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_RuleDResults_EmpiricalVC.csv", header = T, stringsAsFactors = F)
QV_RuleDResults_5050 <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_RuleDResults_FiftyFifty.csv", header = T, stringsAsFactors = F)
QV_RuleEResults_Unif <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_RuleEResults_UniformVC.csv", header = T, stringsAsFactors = F)
QV_RuleEResults_Emp <- read.csv("SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/QV_RuleEResults_EmpiricalVC.csv", header = T, stringsAsFactors = F)

QV_RuleAResults <- mutate(QV_RuleAResults, Rule = "A")
QV_RuleBResults <- mutate(QV_RuleBResults, Rule = "B")
QV_RuleBResults_CNV <- mutate(QV_RuleBResults_CNV, Rule = "B - CNV")
QV_RuleCResults <- mutate(QV_RuleCResults, Rule = "C")
QV_RuleDResults_Unif <- mutate(QV_RuleDResults_Unif, Rule = "D - Unif")
QV_RuleDResults_Emp <- mutate(QV_RuleDResults_Emp, Rule = "D - Emp")
QV_RuleDResults_5050 <- mutate(QV_RuleDResults_5050, Rule = "D - 50/50")
QV_RuleEResults_Unif <- mutate(QV_RuleEResults_Unif, Rule = "E - Unif")
QV_RuleEResults_Emp <- mutate(QV_RuleEResults_Emp, Rule = "E - Emp")

QV_AllRules <- bind_rows(QV_RuleAResults, QV_RuleBResults, QV_RuleBResults_CNV, QV_RuleCResults, QV_RuleDResults_Emp, QV_RuleDResults_Unif, QV_RuleDResults_5050, QV_RuleEResults_Emp, QV_RuleEResults_Unif)




####################################
###### Determine Minority Win ######
###################################%


QV_AllRules <- QV_AllRules %>% 
  mutate(QV1_MinorityWon = ifelse(QV1_Winner_BilingualEduc == "Minority" | QV1_Winner_Immigration == "Minority" | 
                                    QV1_Winner_PublicBonds == "Minority" | QV1_Winner_TeacherTenure == "Minority",T, F),
         QV2_MinorityWon = ifelse(QV2_Winner_BilingualEduc == "Minority" | QV2_Winner_Immigration == "Minority" | 
                                    QV2_Winner_PublicBonds == "Minority" | QV2_Winner_TeacherTenure == "Minority",T, F))


QV_RuleAResults <- filter(QV_AllRules, Rule == "A")
QV_RuleBResults <- filter(QV_AllRules, Rule == "B")
QV_RuleBResults_CNV <- filter(QV_AllRules, Rule == "B - CNV")
QV_RuleCResults <- filter(QV_AllRules, Rule == "C")
QV_RuleDResults_Unif <- filter(QV_AllRules, Rule == "D - Unif")
QV_RuleDResults_Emp <- filter(QV_AllRules, Rule == "D - Emp")
QV_RuleDResults_5050 <- filter(QV_AllRules, Rule == "D - 50/50")
QV_RuleEResults_Unif <- filter(QV_AllRules, Rule == "E - Unif")
QV_RuleEResults_Emp <- filter(QV_AllRules, Rule == "E - Emp")

bootstrapIterations <- nrow(QV_RuleAResults)



###################################################
###### Inequality Measure - Gini Coefficient ######
##################################################%

getGiniCoeff <- function(sampleData, rule, vote, util_NQV, util_QV)
{
  tempRule <- filter(sampleData, Rule == rule)
  
  # get list of simulations in which minority won
  roundsWithMinVic <- tempRule[tempRule[[paste(vote,"_MinorityWon", sep = "")]] == TRUE,]
  roundsWithMinVic <- roundsWithMinVic %>% select(Round) %>% mutate(MinVic = T)
  
  # get utilities under QV and Simple Maj when minority won with QV
  tempRule_UtilNQV <- inner_join(select(tempRule, Round), util_NQV, by = "Round") %>% select(-matches("Round"))
  tempRule_UtilQV <- inner_join(select(tempRule, Round), util_QV, by = "Round") %>% select(-matches("Round"))
  
  # create dataframe to store gini coeff
  tempRule_GiniCoeffs <- data.frame(Round = tempRule$Round, Rule = rule, Gini_UtilNQV = numeric(nrow(tempRule)), Gini_UtilQV = numeric(nrow(tempRule)))
  
  # iterate over each simulation in which minority won to get gini coeff
  for(row in 1:nrow(tempRule))
  {
    print(row)
    tempRule_GiniCoeffs$Gini_UtilNQV[row] <- Gini(tempRule_UtilNQV[row,], corr = FALSE)
    tempRule_GiniCoeffs$Gini_UtilQV[row] <- Gini(tempRule_UtilQV[row,], corr = FALSE)
  }
  
  # Join with RoundsWithMinVic 
  tempRule_GiniCoeffs <- full_join(tempRule_GiniCoeffs, roundsWithMinVic, by = "Round")
  tempRule_GiniCoeffs[is.na(tempRule_GiniCoeffs)] <- F 
  
  # Calculate percentage Gini change between NQV and QV1 & and whether Gini decreased or not
  tempRule_GiniCoeffs <- mutate(tempRule_GiniCoeffs, Gini_UtilNQV = round(Gini_UtilNQV, digits = 4), Gini_UtilQV = round(Gini_UtilQV, digits = 4)) %>%
    mutate(PercentGiniChange = (Gini_UtilQV - Gini_UtilNQV)/Gini_UtilNQV, PercentGiniChange = round(PercentGiniChange, digits = 4),
           ChangeDirection = ifelse(Gini_UtilQV == Gini_UtilNQV, "No Change", ifelse(Gini_UtilQV < Gini_UtilNQV, "Decrease", "Increase")),
           GiniDecrease = ifelse(Gini_UtilQV < Gini_UtilNQV, TRUE, FALSE))
  
  rule <- gsub("/",".", rule)
  
  # Save Gini Raw Data 
  write.csv(tempRule_GiniCoeffs, paste("./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/Ineq - ",vote,"_Rule",rule,"_GiniCoeffs.csv", sep = ""), row.names = F)
  
  return(tempRule_GiniCoeffs)
}

#--- QV1: Calculating Gini Coefficient, percent gini change, & freq of gini decline 

# Calculate Gini coeffs for all sims and all rules

QV1_RuleA_GiniCoeffs <- getGiniCoeff(sampleData = QV_AllRules, rule = "A", vote = "QV1", util_NQV = QV_IndivUtility_NQV_A, util_QV = QV_IndivUtility_QV1_A)
QV1_RuleB_GiniCoeffs <- getGiniCoeff(sampleData = QV_AllRules, rule = "B", vote = "QV1", util_NQV = QV_IndivUtility_NQV_B, util_QV = QV_IndivUtility_QV1_B)
QV1_RuleB_CNV_GiniCoeffs <- getGiniCoeff(sampleData = QV_AllRules, rule = "B - CNV", vote = "QV1", util_NQV = QV_IndivUtility_NQV_B_CNV, util_QV = QV_IndivUtility_QV1_B_CNV)
QV1_RuleC_GiniCoeffs <- getGiniCoeff(sampleData = QV_AllRules, rule = "C", vote = "QV1", util_NQV = QV_IndivUtility_NQV_C, util_QV = QV_IndivUtility_QV1_C)
QV1_RuleD_Unif_GiniCoeffs <- getGiniCoeff(sampleData = QV_AllRules, rule = "D - Unif", vote = "QV1", util_NQV = QV_IndivUtility_NQV_D_Unif, util_QV = QV_IndivUtility_QV1_D_Unif)
QV1_RuleD_Emp_GiniCoeffs <- getGiniCoeff(sampleData = QV_AllRules, rule = "D - Emp", vote = "QV1", util_NQV = QV_IndivUtility_NQV_D_Emp, util_QV = QV_IndivUtility_QV1_D_Emp)
QV1_RuleD_5050_GiniCoeffs <- getGiniCoeff(sampleData = QV_AllRules, rule = "D - 50/50", vote = "QV1", util_NQV = QV_IndivUtility_NQV_D_5050, util_QV = QV_IndivUtility_QV1_D_5050)
QV1_RuleE_Unif_GiniCoeffs <- getGiniCoeff(sampleData = QV_AllRules, rule = "E - Unif", vote = "QV1", util_NQV = QV_IndivUtility_NQV_E_Unif, util_QV = QV_IndivUtility_QV1_E_Unif)
QV1_RuleE_Emp_GiniCoeffs <- getGiniCoeff(sampleData = QV_AllRules, rule = "E - Emp", vote = "QV1", util_NQV = QV_IndivUtility_NQV_E_Emp, util_QV = QV_IndivUtility_QV1_E_Emp)


#--- Modified ---#

# QV1_RuleA_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/Ineq - QV1_RuleA_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# QV1_RuleB_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/Ineq - QV1_RuleB_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# QV1_RuleB_CNV_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/Ineq - QV1_RuleB - CNV_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# QV1_RuleC_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/Ineq - QV1_RuleC_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# QV1_RuleD_Unif_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/Ineq - QV1_RuleD - Unif_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# QV1_RuleD_Emp_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/Ineq - QV1_RuleD - Emp_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# QV1_RuleD_5050_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/Ineq - QV1_RuleD - 50.50_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# QV1_RuleE_Unif_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/Ineq - QV1_RuleE - Unif_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# QV1_RuleE_Emp_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/Ineq - QV1_RuleE - Emp_GiniCoeffs.csv", stringsAsFactors = F, header = T)

meanGiniSummary <- data.frame()

getNthEmpPercentile <- function(nums, NthPercent)
{
  orderedNums <- sort(nums)
  totalNums <- length(orderedNums)
  
  closestPercentile <- round((totalNums/100)*NthPercent)
  
  return(orderedNums[closestPercentile])
}

# Gini of Average Realized Utilities (over all 10k sims)

QV1_AllRules_GiniCoeffs <- bind_rows(QV1_RuleA_GiniCoeffs, QV1_RuleB_GiniCoeffs, QV1_RuleB_CNV_GiniCoeffs, QV1_RuleC_GiniCoeffs, 
                                     QV1_RuleD_Unif_GiniCoeffs, QV1_RuleD_Emp_GiniCoeffs, QV1_RuleD_5050_GiniCoeffs, QV1_RuleE_Unif_GiniCoeffs, QV1_RuleE_Emp_GiniCoeffs)

meanGiniSummary <- bind_rows(meanGiniSummary, QV1_AllRules_GiniCoeffs %>% group_by(Rule) %>% 
                               summarize(MeanGini_Maj = mean(Gini_UtilNQV), lower95CI_Maj = getNthEmpPercentile(Gini_UtilNQV, 2.5), upper95CI_Maj = getNthEmpPercentile(Gini_UtilNQV, 97.5),
                                         MeanGini_QV = mean(Gini_UtilQV), lower95CI_QV = getNthEmpPercentile(Gini_UtilQV, 2.5), upper95CI_QV = getNthEmpPercentile(Gini_UtilQV, 97.5)) %>%
                               mutate(Sample = "QV1", AveragedOver = "All Sims"))


# Filter data so we only make plots conditional on a minority victory

QV1_RuleA_GiniCoeffs <- QV1_RuleA_GiniCoeffs %>% filter(MinVic == T) 
QV1_RuleB_GiniCoeffs <- QV1_RuleB_GiniCoeffs %>% filter(MinVic == T) 
QV1_RuleB_CNV_GiniCoeffs <- QV1_RuleB_CNV_GiniCoeffs %>% filter(MinVic == T) 
QV1_RuleC_GiniCoeffs <- QV1_RuleC_GiniCoeffs %>% filter(MinVic == T) 
QV1_RuleD_Unif_GiniCoeffs <- QV1_RuleD_Unif_GiniCoeffs %>% filter(MinVic == T) 
QV1_RuleD_Emp_GiniCoeffs <- QV1_RuleD_Emp_GiniCoeffs %>% filter(MinVic == T)
QV1_RuleD_5050_GiniCoeffs <- QV1_RuleD_5050_GiniCoeffs %>% filter(MinVic == T) 
QV1_RuleE_Unif_GiniCoeffs <- QV1_RuleE_Unif_GiniCoeffs %>% filter(MinVic == T) 
QV1_RuleE_Emp_GiniCoeffs <- QV1_RuleE_Emp_GiniCoeffs %>% filter(MinVic == T)


# Gini of Average Realized Utilities (over sims with min vic)

QV1_AllRules_GiniCoeffs <- bind_rows(QV1_RuleA_GiniCoeffs, QV1_RuleB_GiniCoeffs, QV1_RuleB_CNV_GiniCoeffs, QV1_RuleC_GiniCoeffs, 
                                     QV1_RuleD_Unif_GiniCoeffs, QV1_RuleD_Emp_GiniCoeffs, QV1_RuleD_5050_GiniCoeffs, QV1_RuleE_Unif_GiniCoeffs, QV1_RuleE_Emp_GiniCoeffs)

meanGiniSummary <- bind_rows(meanGiniSummary, QV1_AllRules_GiniCoeffs %>% group_by(Rule) %>% 
                               summarize(MeanGini_Maj = mean(Gini_UtilNQV), lower95CI_Maj = getNthEmpPercentile(Gini_UtilNQV, 2.5), upper95CI_Maj = getNthEmpPercentile(Gini_UtilNQV, 97.5),
                                         MeanGini_QV = mean(Gini_UtilQV), lower95CI_QV = getNthEmpPercentile(Gini_UtilQV, 2.5), upper95CI_QV = getNthEmpPercentile(Gini_UtilQV, 97.5)) %>%
                               mutate(Sample = "QV1", AveragedOver = "Miv Vic"))


#--- Modified (end) ---#


# QV1 - Summary of Gini Changes

QV1_AllRules_GiniCoeffs <- bind_rows(QV1_RuleA_GiniCoeffs, QV1_RuleB_GiniCoeffs, QV1_RuleB_CNV_GiniCoeffs, QV1_RuleC_GiniCoeffs, QV1_RuleD_Unif_GiniCoeffs, QV1_RuleD_Emp_GiniCoeffs, QV1_RuleD_5050_GiniCoeffs, QV1_RuleE_Unif_GiniCoeffs, QV1_RuleE_Emp_GiniCoeffs)

rulesSummary <- QV1_AllRules_GiniCoeffs %>% group_by(Rule, ChangeDirection) %>% summarize(Count = n()) %>% mutate(Prop = round(Count/sum(Count), digits = 4))

stargazer(rulesSummary, 
          title = "Summary of Gini Changes by Rule (QV1)", type = "text", digits = 3, summary = FALSE, rownames = FALSE, 
          out = "./SV & QV Recalibration on Pop Imm/QV - Inequality/Table - Summary of Gini Changes by Rule (QV1).txt")



# QV1 - Dot Plot of Frequency of Gini Declines per Rule in (Fig 8 - top) 
FreqGiniDeclineA <- group_by(QV1_RuleA_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineB <- group_by(QV1_RuleB_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineB_CNV <- group_by(QV1_RuleB_CNV_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineC <- group_by(QV1_RuleC_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineD_Unif <- group_by(QV1_RuleD_Unif_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineD_Emp <- group_by(QV1_RuleD_Emp_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineD_5050 <- group_by(QV1_RuleD_5050_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineE_Unif <- group_by(QV1_RuleE_Unif_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineE_Emp <- group_by(QV1_RuleE_Emp_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)


QV1_FreqGiniDeclines <- data.frame(Rule = c("A", "B", "B - CNV", "C", "D - Emp", "D - Unif", "D - 50/50","E - Emp", "E - Unif"), 
                                   Freq = c(FreqGiniDeclineA$Freq[1], FreqGiniDeclineB$Freq[1], FreqGiniDeclineB_CNV$Freq[1], 
                                            FreqGiniDeclineC$Freq[1], FreqGiniDeclineD_Emp$Freq[1], FreqGiniDeclineD_Unif$Freq[1], FreqGiniDeclineD_5050$Freq[1],
                                            FreqGiniDeclineE_Emp$Freq[1], FreqGiniDeclineE_Unif$Freq[1]))
QV1_FreqGiniDeclines <- mutate(QV1_FreqGiniDeclines, Freq = round(Freq, digits = 3)*100)

inequalityPlot <- ggplot(data = QV1_FreqGiniDeclines, aes(x=Rule, y = Freq)) + geom_point() + 
  geom_text(aes(label=Freq), hjust = .5, vjust = 2) + ggtitle("QV1 - Percent of Gini Declines by Voting Rule") + 
  scale_y_continuous(limits = c(0,100), breaks = seq(0, 100, by = 10)) +  
  xlab("Bonus Vote - Rules") + ylab("Freq. of Gini Declines (%)") + 
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12))
print(inequalityPlot)
ggsave(inequalityPlot, file="SV & QV Recalibration on Pop Imm/QV - Inequality/QV1 - Percent of Gini Declines by Voting Rule.png", height = 5.5, width = 7.5)



# QV1 - Bar Plot of Mean Percent Gini Change per Rule in (Fig 8 - bottom) 
QV1_MeanGiniCoeffs <- data.frame(Rule = c("A", "B", "B - CNV","C", "D - Emp", "D - Unif", "D - 50/50","E - Emp", "E - Unif"),
                                 MeanPercentGiniChange = c(mean(QV1_RuleA_GiniCoeffs$PercentGiniChange),
                                                           mean(QV1_RuleB_GiniCoeffs$PercentGiniChange),
                                                           mean(QV1_RuleB_CNV_GiniCoeffs$PercentGiniChange),
                                                           mean(QV1_RuleC_GiniCoeffs$PercentGiniChange),
                                                           mean(QV1_RuleD_Emp_GiniCoeffs$PercentGiniChange),
                                                           mean(QV1_RuleD_Unif_GiniCoeffs$PercentGiniChange),
                                                           mean(QV1_RuleD_5050_GiniCoeffs$PercentGiniChange),
                                                           mean(QV1_RuleE_Emp_GiniCoeffs$PercentGiniChange),
                                                           mean(QV1_RuleE_Unif_GiniCoeffs$PercentGiniChange)))

QV1_MeanGiniCoeffs <- mutate(QV1_MeanGiniCoeffs, MeanPercentGiniChange = round(MeanPercentGiniChange, digits = 3)*100)

inequalityPlot <- ggplot(data = QV1_MeanGiniCoeffs, aes(x=Rule, y = MeanPercentGiniChange)) +
  geom_bar(stat="identity") + scale_y_continuous(limits = c(-30,10), breaks = seq(-30, 10, by = 5)) +
  xlab("Bonus Vote - Rules") + ylab("Mean Percent Gini Change") +
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12)) +
  ggtitle("QV1 - Mean Percent Gini Change")
print(inequalityPlot)
ggsave(inequalityPlot, file="SV & QV Recalibration on Pop Imm/QV - Inequality/QV1 - Mean Percent Gini Change.png", height = 5.5, width = 7.5)



# QV1 -  Distribution of Percent Gini Changes

QV1_GiniCoeffs_AllRules <- bind_rows(QV1_RuleA_GiniCoeffs,QV1_RuleB_GiniCoeffs,QV1_RuleB_CNV_GiniCoeffs,QV1_RuleC_GiniCoeffs,QV1_RuleD_Unif_GiniCoeffs,
                                     QV1_RuleD_Emp_GiniCoeffs,QV1_RuleD_5050_GiniCoeffs,QV1_RuleE_Unif_GiniCoeffs,QV1_RuleE_Emp_GiniCoeffs) %>%
  mutate(PercentGiniChange = PercentGiniChange * 100)


inequalityPlot <- ggplot(data = QV1_GiniCoeffs_AllRules, aes(x = PercentGiniChange, colour = Rule)) +
  geom_freqpoly(binwidth = 3) + xlab("Percent Gini Change") + ylab("Absolute Count") +
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12)) +
  ggtitle("QV1 - Smoothed Density of Percent Gini Change (per Rule)")
print(inequalityPlot)
ggsave(inequalityPlot, file="SV & QV Recalibration on Pop Imm/QV - Inequality/QV1 - Smoothed Density of Percent Gini Changes.png", height = 4.5, width = 8.5)







#--- QV2: Calculating Gini Coefficient, percent gini change, & freq of gini decline 

# Calculate Gini coeffs for all sims and all rules

QV2_RuleA_GiniCoeffs <- getGiniCoeff(sampleData = QV_AllRules, rule = "A", vote = "QV2", util_NQV = QV_IndivUtility_NQV_A, util_QV = QV_IndivUtility_QV2_A)
QV2_RuleB_GiniCoeffs <- getGiniCoeff(sampleData = QV_AllRules, rule = "B", vote = "QV2", util_NQV = QV_IndivUtility_NQV_B, util_QV = QV_IndivUtility_QV2_B)
QV2_RuleB_CNV_GiniCoeffs <- getGiniCoeff(sampleData = QV_AllRules, rule = "B - CNV", vote = "QV2", util_NQV = QV_IndivUtility_NQV_B_CNV, util_QV = QV_IndivUtility_QV2_B_CNV)
QV2_RuleC_GiniCoeffs <- getGiniCoeff(sampleData = QV_AllRules, rule = "C", vote = "QV2", util_NQV = QV_IndivUtility_NQV_C, util_QV = QV_IndivUtility_QV2_C)
QV2_RuleD_Unif_GiniCoeffs <- getGiniCoeff(sampleData = QV_AllRules, rule = "D - Unif", vote = "QV2", util_NQV = QV_IndivUtility_NQV_D_Unif, util_QV = QV_IndivUtility_QV2_D_Unif)
QV2_RuleD_Emp_GiniCoeffs <- getGiniCoeff(sampleData = QV_AllRules, rule = "D - Emp", vote = "QV2", util_NQV = QV_IndivUtility_NQV_D_Emp, util_QV = QV_IndivUtility_QV2_D_Emp)
QV2_RuleD_5050_GiniCoeffs <- getGiniCoeff(sampleData = QV_AllRules, rule = "D - 50/50", vote = "QV2", util_NQV = QV_IndivUtility_NQV_D_5050, util_QV = QV_IndivUtility_QV2_D_5050)
QV2_RuleE_Unif_GiniCoeffs <- getGiniCoeff(sampleData = QV_AllRules, rule = "E - Unif", vote = "QV2", util_NQV = QV_IndivUtility_NQV_E_Unif, util_QV = QV_IndivUtility_QV2_E_Unif)
QV2_RuleE_Emp_GiniCoeffs <- getGiniCoeff(sampleData = QV_AllRules, rule = "E - Emp", vote = "QV2", util_NQV = QV_IndivUtility_NQV_E_Emp, util_QV = QV_IndivUtility_QV2_E_Emp)


#--- Modified ---#

# QV2_RuleA_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/Ineq - QV2_RuleA_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# QV2_RuleB_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/Ineq - QV2_RuleB_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# QV2_RuleB_CNV_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/Ineq - QV2_RuleB - CNV_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# QV2_RuleC_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/Ineq - QV2_RuleC_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# QV2_RuleD_Unif_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/Ineq - QV2_RuleD - Unif_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# QV2_RuleD_Emp_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/Ineq - QV2_RuleD - Emp_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# QV2_RuleD_5050_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/Ineq - QV2_RuleD - 50.50_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# QV2_RuleE_Unif_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/Ineq - QV2_RuleE - Unif_GiniCoeffs.csv", stringsAsFactors = F, header = T)
# QV2_RuleE_Emp_GiniCoeffs <- read.csv("./SV & QV Recalibration on Pop Imm/QV - Raw Simulation Data/Ineq - QV2_RuleE - Emp_GiniCoeffs.csv", stringsAsFactors = F, header = T)


# Gini of Average Realized Utilities (over all 10k sims)

QV2_AllRules_GiniCoeffs <- bind_rows(QV2_RuleA_GiniCoeffs, QV2_RuleB_GiniCoeffs, QV2_RuleB_CNV_GiniCoeffs, QV2_RuleC_GiniCoeffs, QV2_RuleD_Unif_GiniCoeffs, QV2_RuleD_Emp_GiniCoeffs, QV2_RuleD_5050_GiniCoeffs, QV2_RuleE_Unif_GiniCoeffs, QV2_RuleE_Emp_GiniCoeffs)

meanGiniSummary <- bind_rows(meanGiniSummary, QV2_AllRules_GiniCoeffs %>% group_by(Rule) %>% 
                               summarize(MeanGini_Maj = mean(Gini_UtilNQV), lower95CI_Maj = getNthEmpPercentile(Gini_UtilNQV, 2.5), upper95CI_Maj = getNthEmpPercentile(Gini_UtilNQV, 97.5),
                                         MeanGini_QV = mean(Gini_UtilQV), lower95CI_QV = getNthEmpPercentile(Gini_UtilQV, 2.5), upper95CI_QV = getNthEmpPercentile(Gini_UtilQV, 97.5)) %>%
                               mutate(Sample = "QV2", AveragedOver = "All Sims"))


# Filter data so we only make plots conditional on a minority victory

QV2_RuleA_GiniCoeffs <- QV2_RuleA_GiniCoeffs %>% filter(MinVic == T) 
QV2_RuleB_GiniCoeffs <- QV2_RuleB_GiniCoeffs %>% filter(MinVic == T) 
QV2_RuleB_CNV_GiniCoeffs <- QV2_RuleB_CNV_GiniCoeffs %>% filter(MinVic == T) 
QV2_RuleC_GiniCoeffs <- QV2_RuleC_GiniCoeffs %>% filter(MinVic == T) 
QV2_RuleD_Unif_GiniCoeffs <- QV2_RuleD_Unif_GiniCoeffs %>% filter(MinVic == T) 
QV2_RuleD_Emp_GiniCoeffs <- QV2_RuleD_Emp_GiniCoeffs %>% filter(MinVic == T)
QV2_RuleD_5050_GiniCoeffs <- QV2_RuleD_5050_GiniCoeffs %>% filter(MinVic == T) 
QV2_RuleE_Unif_GiniCoeffs <- QV2_RuleE_Unif_GiniCoeffs %>% filter(MinVic == T) 
QV2_RuleE_Emp_GiniCoeffs <- QV2_RuleE_Emp_GiniCoeffs %>% filter(MinVic == T)


# Gini of Average Realized Utilities (over sims with min vic)

QV2_AllRules_GiniCoeffs <- bind_rows(QV2_RuleA_GiniCoeffs, QV2_RuleB_GiniCoeffs, QV2_RuleB_CNV_GiniCoeffs, QV2_RuleC_GiniCoeffs, QV2_RuleD_Unif_GiniCoeffs, QV2_RuleD_Emp_GiniCoeffs, QV2_RuleD_5050_GiniCoeffs, QV2_RuleE_Unif_GiniCoeffs, QV2_RuleE_Emp_GiniCoeffs)

meanGiniSummary <- bind_rows(meanGiniSummary, QV2_AllRules_GiniCoeffs %>% group_by(Rule) %>% 
                               summarize(MeanGini_Maj = mean(Gini_UtilNQV), lower95CI_Maj = getNthEmpPercentile(Gini_UtilNQV, 2.5), upper95CI_Maj = getNthEmpPercentile(Gini_UtilNQV, 97.5),
                                         MeanGini_QV = mean(Gini_UtilQV), lower95CI_QV = getNthEmpPercentile(Gini_UtilQV, 2.5), upper95CI_QV = getNthEmpPercentile(Gini_UtilQV, 97.5)) %>%
                               mutate(Sample = "QV2", AveragedOver = "Miv Vic"))

write.csv(meanGiniSummary, "./SV & QV Recalibration on Pop Imm/QV - Inequality/Mean Gini Coeffs & CIs.csv", row.names = F)

#--- Modified (end) ---#



# Filter data so we only make plots conditional on a minority victory


# QV2 - Summary of Gini Changes

QV2_AllRules_GiniCoeffs <- bind_rows(QV2_RuleA_GiniCoeffs, QV2_RuleB_GiniCoeffs, QV2_RuleB_CNV_GiniCoeffs, QV2_RuleC_GiniCoeffs, QV2_RuleD_Unif_GiniCoeffs, QV2_RuleD_Emp_GiniCoeffs, QV2_RuleD_5050_GiniCoeffs, QV2_RuleE_Unif_GiniCoeffs, QV2_RuleE_Emp_GiniCoeffs)

rulesSummary <- QV2_AllRules_GiniCoeffs %>% group_by(Rule, ChangeDirection) %>% summarize(Count = n()) %>% mutate(Prop = round(Count/sum(Count), digits = 4))

stargazer(rulesSummary, 
          title = "Summary of Gini Changes by Rule (QV2)", type = "text", digits = 3, summary = FALSE, rownames = FALSE, 
          out = "./SV & QV Recalibration on Pop Imm/QV - Inequality/Table - Summary of Gini Changes by Rule (QV2).txt")



# QV2 - Dot Plot of Frequency of Gini Declines per Rule in (Fig 8 - top) 
FreqGiniDeclineA <- group_by(QV2_RuleA_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineB <- group_by(QV2_RuleB_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineB_CNV <- group_by(QV2_RuleB_CNV_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineC <- group_by(QV2_RuleC_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineD_Unif <- group_by(QV2_RuleD_Unif_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineD_Emp <- group_by(QV2_RuleD_Emp_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineD_5050 <- group_by(QV2_RuleD_5050_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineE_Unif <- group_by(QV2_RuleE_Unif_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)
FreqGiniDeclineE_Emp <- group_by(QV2_RuleE_Emp_GiniCoeffs, GiniDecrease) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) %>% filter(GiniDecrease == T)


QV2_FreqGiniDeclines <- data.frame(Rule = c("A", "B", "B - CNV","C", "D - Emp", "D - Unif", "D - 50/50","E - Emp", "E - Unif"), 
                                   Freq = c(FreqGiniDeclineA$Freq[1], FreqGiniDeclineB$Freq[1], FreqGiniDeclineB_CNV$Freq[1], FreqGiniDeclineC$Freq[1], 
                                            FreqGiniDeclineD_Emp$Freq[1], FreqGiniDeclineD_Unif$Freq[1], FreqGiniDeclineD_5050$Freq[1],
                                            FreqGiniDeclineE_Emp$Freq[1], FreqGiniDeclineE_Unif$Freq[1]))
QV2_FreqGiniDeclines <- mutate(QV2_FreqGiniDeclines, Freq = round(Freq, digits = 3)*100)

inequalityPlot <- ggplot(data = QV2_FreqGiniDeclines, aes(x=Rule, y = Freq)) + geom_point() + 
  geom_text(aes(label=Freq), hjust = .5, vjust = 2) + ggtitle("QV2 - Percent of Gini Declines by Voting Rule") + 
  scale_y_continuous(limits = c(0,100), breaks = seq(0, 100, by = 10)) +  
  xlab("Bonus Vote - Rules") + ylab("Freq. of Gini Declines (%)") + 
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12))
print(inequalityPlot)
ggsave(inequalityPlot, file="SV & QV Recalibration on Pop Imm/QV - Inequality/QV2 - Percent of Gini Declines by Voting Rule.png", height = 5.5, width = 7.5)



# QV2 - Bar Plot of Mean Percent Gini Change per Rule in (Fig 8 - bottom) 
QV2_MeanGiniCoeffs <- data.frame(Rule = c("A", "B", "B - CNV", "C", "D - Emp", "D - Unif", "D - 50/50","E - Emp", "E - Unif"),
                                 MeanPercentGiniChange = c(mean(QV2_RuleA_GiniCoeffs$PercentGiniChange),
                                                           mean(QV2_RuleB_GiniCoeffs$PercentGiniChange),
                                                           mean(QV2_RuleB_CNV_GiniCoeffs$PercentGiniChange),
                                                           mean(QV2_RuleC_GiniCoeffs$PercentGiniChange),
                                                           mean(QV2_RuleD_Emp_GiniCoeffs$PercentGiniChange),
                                                           mean(QV2_RuleD_Unif_GiniCoeffs$PercentGiniChange),
                                                           mean(QV2_RuleD_5050_GiniCoeffs$PercentGiniChange),
                                                           mean(QV2_RuleE_Emp_GiniCoeffs$PercentGiniChange),
                                                           mean(QV2_RuleE_Unif_GiniCoeffs$PercentGiniChange)))

QV2_MeanGiniCoeffs <- mutate(QV2_MeanGiniCoeffs, MeanPercentGiniChange = round(MeanPercentGiniChange, digits = 3)*100)

inequalityPlot <- ggplot(data = QV2_MeanGiniCoeffs, aes(x=Rule, y = MeanPercentGiniChange)) +
  geom_bar(stat="identity") + scale_y_continuous(limits = c(-30,10), breaks = seq(-30, 10, by = 5)) +
  xlab("Bonus Vote - Rules") + ylab("Mean Percent Gini Change") +
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12)) +
  ggtitle("QV2 - Mean Percent Gini Change")
print(inequalityPlot)
ggsave(inequalityPlot, file="SV & QV Recalibration on Pop Imm/QV - Inequality/QV2 - Mean Percent Gini Change.png", height = 5.5, width = 7.5)



# QV2 -  Distribution of Percent Gini Changes

QV2_GiniCoeffs_AllRules <- bind_rows(QV2_RuleA_GiniCoeffs,QV2_RuleB_GiniCoeffs,QV2_RuleB_CNV_GiniCoeffs,QV2_RuleC_GiniCoeffs,QV2_RuleD_Unif_GiniCoeffs,
                                     QV2_RuleD_Emp_GiniCoeffs,QV2_RuleD_5050_GiniCoeffs,QV2_RuleE_Unif_GiniCoeffs,QV2_RuleE_Emp_GiniCoeffs) %>%
  mutate(PercentGiniChange = PercentGiniChange * 100)


inequalityPlot <- ggplot(data = QV2_GiniCoeffs_AllRules, aes(x = PercentGiniChange, colour = Rule)) +
  geom_freqpoly(binwidth = 3) + xlab("Percent Gini Change") + ylab("Absolute Count") +
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12)) +
  ggtitle("QV2 - Smoothed Density of Percent Gini Change (per Rule)")
print(inequalityPlot)
ggsave(inequalityPlot, file="SV & QV Recalibration on Pop Imm/QV - Inequality/QV2 - Smoothed Density of Percent Gini Changes.png", height = 4.5, width = 8.5)



